TP 4 : Stratification

Auteur·rice

Paul Géhin

Date de publication

27 mars 2026

Exercice 1

Le dirigeant d’un pays de 100 000 habitants souhaite connaître la fréquentation des théâtres au cours des 12 derniers mois. Pour faire cela en respectant des contraintes de coût, il fait appel à un statisticien lui recommandant d’estimer cet indicateur à l’aide d’un sondage. Le professionnel dispose d’une base de sondage des individus du pays dans laquelle se trouve deux informations auxiliaires :

  • dep : le département dans lequel l’individu vit.
  • CSP : la catégorie socio-professionnelle de l’individu en 10 classes (il ne s’agit pas de la nomenclature PCS 2020 de l’Insee).

On suppose que les théâtres sont répartis de manière similaire dans tous les départements et que tous les individus peuvent y accéder dans les mêmes conditions.

Compte tenu des coûts, l’échantillon ne peut contenir plus de 10 000 individus.

Le statisticien décide alors de comparer trois approches :

  • tirer un échantillon de 10 000 individus à l’aide d’un SRS parmi l’ensemble des habitants du pays. (tirage 1)
  • tirer un échantillon de 10 000 individus à l’aide d’un SRS stratifié par département avec allocation proportionnelle.(tirage 2)
  • tirer un échantillon de 10 000 individus à l’aide d’un SRSR stratifié par CSP avec allocation proportionnelle. (tirage 3)

Les résultats des tirages sont stockés dans les tables ech1.csv (pour le tirage 1), ech2.csv (pour le tirage 2) et ech3 (pour le tirage 3) et contiennent les informations auxiliaires ainsi qu’une variable theatre indiquant le nombre de visites au théâtre au cours des 12 derniers mois et une variable Prob indiquant la probabilité d’inclusion d’ordre un.

Dans la suite, on s’intéresse à la fréquentation moyenne par habitant et on note \(y_k\) le nombre de fois où l’individu \(k\) est allé au théâtre au cours des 12 derniers mois.

  1. Décrivez la population, la variable d’intérêt et la fonction d’intérêt. Dans la suite, nous noterons, \(\mathcal{U}\) , la population.
  • La population est l’ensemble des individus du pays.
  • La variable d’intérêt correspond au nombre de fréquentation aux théâtres au cours des 12 derniers mois.
  • La fonction d’intérêt est la moyenne des fréquentations.
  1. Donnez un estimateur sans biais de \(\displaystyle \mu_y = \frac{1}{N} \sum_{k \in \mathcal{U}} y_k\) basé sur l’estimateur d’Horvitz-Thompson et la taille de la population \(N\).

Un estimateur sans biais de \(\mu_y\) noté \(\hat{\mu}_y\) est donné par \(\displaystyle \hat{\mu}_y = \frac{1}{N} \hat{t}_{y,\text{HT}}\). En effet, en utilisant la linéarité de l’opérateur d’espérance, il en vient que \(\hat{\mu}_y\) est un estimateur sans biais de \(\mu_y\)

Cas du tirage 1

  1. Proposez une estimation associée à cet estimateur.
library("sampling")

ech1 <- read.csv("https://sondages.cours.gehin.net/TP/TP%204/data/ech1.csv", colClasses = c("dep" = "character"))
n <- nrow(ech1)
N <- 100000
#Vérification du contenu de la table
str(ech1)

#Moyenne
#Il suffit de diviser notre estimation du total par la taille de la population
estim_moy_1 <- HTestimator(ech1$theatre, ech1$Prob)/N

cat("L'estimation de la moyenne associée est :", estim_moy_1)
library("sampling")

ech1 <- read.csv("https://sondages.cours.gehin.net/TP/TP%204/data/ech1.csv", colClasses = c("dep" = "character"))
n <- nrow(ech1)
N <- 100000
#Vérification du contenu de la table
str(ech1)

#Total
estim1 <- HTestimator(ech1$theatre, ech1$Prob)

cat("L'estimation associée à l'estimateur du total est : ", estim1)
  1. Proposez un estimateur de la variance de l’estimateur proposé à la question 2. On pourra se rappeler que pour une variable aléatoire \(X\) et un réel \(\lambda\), \(\mathbb{V}(\lambda X) = \lambda^2 \mathbb{V}(X)\).

Nous avons \(\displaystyle \hat{\mu}_y = \frac{1}{N} \hat{t}_{y,\text{HT}}\). Il en vient que \(\displaystyle \mathbb{V}(\hat{\mu}_y) = \mathbb{V}(\frac{1}{N} \hat{t}_{y,\text{HT}}) = \frac{1}{N^2} \mathbb{V}(\hat{t}_{y,\text{HT}})\).

Dans le cadre d’un SRS, il est possible d’utiliser l’estimateur HT de la variance ou l’estimateur de Sen-Yates-Grundy indifféremment car ils sont égaux. De plus, sous un SRS, les probabilités d’inclusion d’ordre 2 sont non nulles pour tous les individus de la population \(\mathcal{U}\).

Il en vient que la variance \(\displaystyle \mathbb{V}(\hat{\mu}_y)\) peut être estimée par \(\displaystyle \hat{\mathbb{V}}(\hat{\mu}_y) = \frac{1}{N} \hat{\mathbb{V}}_\text{SYG}(\hat{t}_{y,\text{HT}})\).

  1. Donnez un intervalle de confiance asymptotique au niveau 0.90 du total. Calculez une réalisation de cet intervalle.

Pour calculer un intervalle de confiance asymptotique, il faut : - un estimateur du paramètre. - un estimateur de la variance du paramètre. - les quantiles de la loi normale.

#Calcul d'un estimateur de la variance du total.
#Ici l'estimateur de SYG et HT donnent les mêmes résultats car SRS.

var_moy_1 <- varHT(ech1$theatre, pikl = matrice_prob_double(n,N))/N^2
alpha <- 0.1
q <- qnorm(1- (alpha/2))

#Calcul des bornes de la réalisation de l'invervalle.
binf_moy_1 <- estim_moy_1 - q*sqrt(var_moy_1)
bsup_moy_1 <- estim_moy_1 + q*sqrt(var_moy_1)


cat("Une réalisation de l'invervalle de confiance au niveau 1-alpha est [", binf_moy_1, " ; ", bsup_moy_1, "]")
#Calcul d'un estimateur de la variance du total.
#Ici l'estimateur de SYG et HT donnent les mêmes résultats car SRS.

var1 <- varHT(ech1$theatre, pikl = matrice_prob_double(n,N))
alpha <- 0.1
q <- qnorm(1- (alpha/2))

#Calcul des bornes de la réalisation de l'invervalle.
binf1 <- estim1 - q*sqrt(var1)
bsup1 <- estim1 + q*sqrt(var1)

cat("Une réalisation de l'invervalle de confiance au niveau 1-alpha est [", binf1, " ; ", bsup1, "]")

Cas du tirage 2

  1. En remarquant que \(\hat{t}_{y,\text{HT}} = \sum_{h \in [H]} \hat{t}_{y,\text{HT},h}\) où \([H] = \{1, ..., H\}\) désigne l’ensemble des strates associées à une stratification et \(\hat{t}_{y,\text{HT},h} = \sum_{k \in \mathcal{s} \cap \mathcal{U}_h} \frac{y_k}{\pi_k}\), réécrivez l’estimateur proposé à la question 2 comme une somme d’estimateurs calculé par département.

Dans cette question (comme dans le cours), nous notons \(h\) les strates. Ici les strates correspondent aux départements du pays.

On sait que \(\hat{\mu}_y = \frac{1}{N} \hat{t}_{y,\text{HT}}\). Or d’après en utilisant la remarque : \(\hat{t}_{y,\text{HT}} = \sum_{h \in [H]} \hat{t}_{y,\text{HT},h}\), on en déduit que \[\begin{align} \hat{\mu}_y &= \frac{1}{N} \hat{t}_{y,\text{HT}} \\ &= \frac{1}{N} \sum_{h \in [H]} \hat{t}_{y,\text{HT},h} \\ &= \sum_{h \in [H]} \frac{1}{N} \hat{t}_{y,\text{HT},h} \\ \end{align}\]

Pour calculer une estimation de la moyenne sur l’ensemble de la population, il suffira de calculer les totaux dans chaque strate et de diviser par la taille de population totale (et non par strate). Nous allons appliquer ce calcul dans la question suivante.

  1. Proposez une estimation associée à cet estimateur.

Le code ci-dessous n’est pas optimisé (utilisation de boucles en R, …) : il vise à être accessible pour des nouveaux utilisateurs dans ce langage.

ech2 <- read.csv("https://sondages.cours.gehin.net/TP/TP%204/data/ech2.csv", colClasses = c("dep" = "character"))
n <- nrow(ech2)
N <- 100000
#Vérification du contenu de la table
str(ech2)


#Récupération de la liste des établissements = liste des strates
liste_dep <- unique(ech2$dep)


#Pour chaque département on calcule le terme dans la somme de la question précédente

estim_moy_2 <- 0

for(i in liste_dep){
    estim_moy_2 <- estim_moy_2 + HTestimator(ech2[ech2$dep == i, "theatre"],
                             ech2[ech2$dep == i, "Prob"])/N
}


cat("L'estimation de la moyenne associée est :", estim_moy_2)
ech2 <- read.csv("Donnees/ech2.csv", colClasses = c("dep" = "character"))
n <- nrow(ech2)
N <- 100000
#Vérification du contenu de la table
str(ech2)

#La fonction HTstrata permet de calculer directement en indiquant les strates
#l'estimation du total. Il suffit de diviser par la taille de la population afin
#d'obtenir un estimateur de la moyenne sur l'ensemble de la population.

estim_tot_2_direct <- HTstrata(y = ech2$theatre, pik = ech2$Prob, strata = ech2$dep)
estim_moy_2_direct <- estim_tot_2_direct/N

cat("L'estimation de la moyenne associée est :", estim_moy_2_direct)
ech2 <- read.csv("https://sondages.cours.gehin.net/TP/TP%204/data/ech2.csv", colClasses = c("dep" = "character"))
n <- nrow(ech2)
N <- 100000
#Vérification du contenu de la table
str(ech2)


#Récupération de la liste des établissements = liste des strates
liste_dep <- unique(ech2$dep)


#Pour chaque département on calcule le total dans chaque strate (= département)
estim_tot_2 <- 0

for(i in liste_dep){
    estim_tot_2 <- estim_tot_2 + HTestimator(ech2[ech2$dep == i, "theatre"],
                             ech2[ech2$dep == i, "Prob"])
}


cat("L'estimation de la moyenne associée est :", estim_tot_2)
  1. Proposez un estimateur de la variance de l’estimateur proposé à la question 6.

La variance d’un estimateur stratifié est égale à la somme des variances des estimateurs resteints à une strate en vertu des tirages indépendants entre strates. Il en vient que la variance de \(\hat{\mu}_2\), \[ \begin{align} \mathbb{V}(\hat{\mu}_y) &= \mathbb{V}( \sum_{h \in [H]} \frac{1}{N} \hat{t}_{y,\text{HT},h}) \text{ en utilisant la question 6} \\ &= \sum_{h \in [H]} \mathbb{V}( \frac{1}{N} \hat{t}_{y,\text{HT},h}) \text{ par indépendance des tirages entre strates} \\ &= \sum_{h \in [H]} \frac{1}{N^2} \mathbb{V}(\hat{t}_{y,\text{HT},h}) \end{align} \]

La variance de l’estimateur stratifié de la moyenne sera donc égal à la somme des variances des estimateurs totaux divisée par la taille de l’échantillon au carré.

Dans le cadre de cet exercice, les tirages au sein des strates est un sondage aléatoire simple sans remise.

La variance de \(\hat{\mu}_2\) notée \(\mathbb{V}(\hat{\mu}_2)\) est donnée par \[ \begin{align} \mathbb{V}(\hat{\mu}_y) &= \sum_{h \in [H]} \frac{1}{N^2} \mathbb{V}( \hat{t}_{y,\text{HT},h}) \\ &= \frac{1}{N^2} \sum_{h \in [H]} \mathbb{V}(\hat{t}_{y,\text{HT},h}) \\ &= \frac{1}{N^2} \sum_{h \in [H]} \frac{N_h^2}{n_h} \left(1 - \frac{n_h}{N_h} \right) S_{y,h}^2 \\ & \text{avec } S_{y,h}^2 = \frac{1}{N_h - 1} \sum_{j \in \mathcal{U}_h} \left(y_k - \bar{y}_h \right)^2 \end{align} \]

Un estimateur sans biais de \(\mathbb{V}(\hat{\mu}_y)\) noté \(\hat{\mathbb{V}}(\hat{\mu}_y)\) (car les probabilités d’inclusion d’ordre 2 sont non nulles pour tout individu de la population) est donné par \[ \begin{align} \hat{\mathbb{V}}(\hat{\mu}_y) &= \frac{1}{N^2} \sum_{h \in [H]} \frac{N_h^2}{n_h} \left(1 - \frac{n_h}{N_h} \right) \hat{S}_{y,h}^2 \\ & \text{avec } \hat{S}_{y,h}^2 = \frac{1}{N_h - 1} \sum_{j \in {S}_h} \left(y_k - \tilde{y}_h \right)^2 \text{ et } \tilde{y}_h = \sum_{k \in S_h} y_k \end{align} \]

  1. Donnez un intervalle de confiance asymptotique au niveau 0.90 du total. Calculez une réalisation de cet intervalle.

On reprend les mêmes étapes qu’habituellement :

  • calcul de l’estimation de la fonction d’intérêt.
  • calcul de l’estimation de la variance de l’estimateur de la fonction d’intérêt.
  • calcul des quantiles de la loi normale.

Pour calculer l’estimation de la variance de l’estimateur, nous avons besoin en vertu de la question précédente de la taille de la population \(N_h\). Or cette information n’est pas directement disponible dans les tables. Cependant, nous disposons de la probabilité d’inclusion d’ordre un pour chaque individu ainsi que de la taille de l’échantillon tiré dans chaque strate (en calculant le nombre d’individus par strate à l’aide de la table ech2).

Or on sait que pour tout individu \(k\) de la strate \(h\), \(\pi_k = \frac{n_h}{N_h}\) donc \(N_h = \pi_k \times n_h\).

#Calcul de n, N par strate et stockage dans allocation_dep
#car nécessaire au calcul de la variance comme indiqué dans la question précédente

n_ech2 <- table(ech2$dep)
prob_ech2 <- unique(ech2[, c("dep", "Prob")])
prob_ech2 <- setNames(prob_ech2$Prob, prob_ech2$dep)

N_ech2 <- n_ech2/prob_ech2

allocation_dep <- data.frame("n" = as.vector(n_ech2),
                             "N" = as.vector(N_ech2),
                             "dep" = names(n_ech2))

var2 <- 0

for(i in allocation_dep$dep){
  var2 <- var2 + sampling::varHT(ech2[ech2$dep == i, "theatre"],
                                 matrice_prob_double(allocation_dep[allocation_dep$dep == i, "n"],
                                                     allocation_dep[allocation_dep$dep == i, "N"]))/N^2
}

binf2 <- estim_moy_2 - q*sqrt(var2)
bsup2 <- estim_moy_2 + q*sqrt(var2)


cat("Un réalisation d'un intervalle de confiance à 90% est donné par : [", binf2*N, ";", bsup2*N, "]")

Version total

Bientôt disponible

:::

  1. Comparez avec le résultat de la question 5. Est-ce que la stratification par département est efficace ? Pourquoi ?

Cas du tirage 3

  1. Donnez l’allocation associée à ce plan stratifié.

L’allocation \((n_1, ..., n_H)\) correspond au nombre d’individus tirés par strate. Ici les strates correspondent aux CSP. Il suffit donc de compter le nombre d’individus dans l’échantillon par CSP.

ech3 <- read.csv("https://sondages.cours.gehin.net/TP/TP%204/data/ech2.csv", colClasses = c("dep" = "character"))
n <- nrow(ech3)
N <- 100000
#Vérification du contenu de la table
str(ech3)

#Calcul de l'allocation
table(ech3$CSP)
  1. Utilisez l’estimateur proposé à la question 6 avec cette nouvelle stratification.

Même réponse qu’à la question 8. en utilisant les CSP à la place des départements.

  1. Proposez une estimation associée à cet estimateur.
ech3 <- read.csv("https://sondages.cours.gehin.net/TP/TP%204/data/ech3.csv", colClasses = c("dep" = "character"))
n <- nrow(ech3)
N <- 100000
#Vérification du contenu de la table
str(ech3)


#Récupération de la liste des établissements = liste des strates
liste_CSP <- unique(ech3$CSP)


#Pour chaque département on calcule le terme dans la somme de la question précédente

estim_moy_3 <- 0

for(i in liste_CSP){
    estim_moy_3 <- estim_moy_3 + HTestimator(ech3[ech3$CSP == i, "theatre"],
                             ech3[ech3$CSP == i, "Prob"])/N
}


cat("L'estimation de la moyenne associée est :", estim_moy_3)
  1. Proposez un estimateur de la variance de l’estimateur proposé à la question 12.

Même réponse qu’à la question 9. en utilisant les CSP à la place des départements.

  1. Donnez un intervalle de confiance asymptotique au niveau 0.90 du total. Calculez une réalisation de cet intervalle.
#Calcul de n, N par strate et stockage dans allocation_dep
#car nécessaire au calcul de la variance comme indiqué dans la question précédente

n_ech3 <- table(ech3$CSP)
prob_ech3 <- unique(ech3[, c("CSP", "Prob")])
prob_ech3 <- setNames(prob_ech3$Prob, prob_ech3$CSP)

N_ech3 <- n_ech3/prob_ech3

allocation_dep <- data.frame("n" = as.vector(n_ech3),
                             "N" = as.vector(N_ech3),
                             "CSP" = names(n_ech3))

var3 <- 0

for(i in allocation_dep$CSP){
  var3 <- var3 + sampling::varHT(ech3[ech3$CSP == i, "theatre"],
                                 matrice_prob_double(allocation_dep[allocation_dep$CSP == i, "n"],
                                                     allocation_dep[allocation_dep$CSP == i, "N"]))/N^2
}

binf3 <- estim_moy_3 - q*sqrt(var3)
bsup3 <- estim_moy_3 + q*sqrt(var3)

cat("Un réalisation d'un intervalle de confiance à 90% est donné par : [", binf3*N, ";", bsup3*N, "]")
  1. Comparez avec le résultat de la question 5 et 10. Est-ce que la stratification par département est efficace ? Pourquoi ?
  • La variance de l’estimateur de moyenne sans stratification vaut r var1/N^2.
  • La variance de l’estimateur de moyenne dans le cas d’une stratification par département vaut r var2.
  • La variance de l’estimateur de moyenne dans le cas d’une stratification par CSP vaut r var3.

La variance de l’estimateur sans stratification et avec stratification par département sont environ les mêmes : ceci tient au fait que le département dans lequel vit l’individu est indépendamment de la fréquence en théâtre (dans cet exercice que le nombre de théâtre et l’accès aux théâtres étaient les mêmes d’un département à un autre.)

Par contre, la variance de l’estimateur par stratification par CSP est en deçà : la CSP explique une partie de la variabilité en fréquentation en théâtre et permet des gains lorsque cette variable est mobilisée pour la stratification.

Attention : ici, nous avons des gains en utilisant la variable CSP car il y a un fort lien entre CSP et fréquentation en théâtre. Pour une autre variable d’intérêt, une stratification par CSP peut être inutile : par exemple, si la variable d’intérêt est indépendante de la CSP.